function [A2, A3] = A_23(NmE, NmF1, A1, hmF2, hmF1, hmE, BEtop, B1bot, B2bot, foF1)
%
% F1-Layer Amplitude (A2) and E-Layer Amplitude (A3)
%
%DESCRIPTION:
%This function computes the F1-Layer Amplitude (A2) and E-Layer
%Amplitude (A3).
%
%PROTOTYPE:
% [A2, A3] = A_23(NmE, NmF1, A1, hmF2, hmF1, hmE, BEtop, B1bot, B2bot, foF1)
%
%--------------------------------------------------------------------------
% INPUTS:
%   NmE        [1x1]       E-Layer Max. Density      [1e11 m-3]
%   NmF1       [1x1]       F1-Layer Max. Density     [1e11 m-3]
%   A1         [1x1]       F2-Layer Amplitude        [1e11 m-3]
%   hmF2       [1x1]       F2-Layer Max. Den. Height [km]
%   hmF1       [1x1]       F1-Layer Max. Den. Height [km]
%   hmE        [1x1]       E-Layer Max. Den. Height  [km]
%   BEtop      [1x1]       E-Layer Top Thickness     [km]
%   B1bot      [1x1]       F1-Layer Bottom Thickness [km]
%   B2bot      [1x1]       F2-Layer Bottom Thickness [km]
%   foF1       [1x1]       F2-Layer Crit. Freq.      [MHz]
%--------------------------------------------------------------------------
% OUTPUTS:
%   A2         [1x1]       F1-Layer Amplitude        [1e11 m-3]
%   A3         [1x1]       E-Layer Amplitude         [1e11 m-3]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

if foF1 < 0.5
    aux = Epst(A1, hmF2, B2bot, hmE);
    A2 = 0;
    A3 = 4.0 * (NmE - aux);
elseif foF1 >= 0.5
    A3a = 4 * NmE;
    for k = 1 : 5
        aux1 = Epst(A1, hmF2, B2bot, hmF1);
        aux2 = Epst(A3a, hmE, BEtop, hmF1);
        A2a = 4 * (NmF1 - aux1 - aux2);
        A2a = (A2a * exp(A2a - 0.8*NmF1) + 0.8*NmF1) / (1 + exp(A2a - 0.8*NmF1));
        aux3 = Epst(A2a, hmF1, B1bot, hmE);
        aux4 = Epst(A1, hmF2, B2bot, hmE);
        A3a = 4 * (NmE - aux3 - aux4);
    end
    A2 = A2a;
    A3 = (A3a*exp(60*(A3a - 0.005)) + 0.05) / (1 + exp(60*(A3a - 0.005)));
end

end